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_^ Abstract 

^ . We describe the new version 2.00c of the code hfbtho that solves the nuclear Skyrme 

! Hartree-Fock (HF) or Skyrme Hartree-Fock-Bogolyubov (HFB) problem by using the cylindri- 
cal transformed deformed harmonic-oscillator basis. In the new version, we have implemented 
the following features: (i) the modified Broyden method for non-linear problems, (ii) optional 
breaking of reflection symmetry, (iii) calculation of axial multipole moments, (iv) finite temper- 
ed ■ ature formalism for the HFB method, (v) linear constraint method based on the approximation 
of the Random Phase Approximation (RPA) matrix for multi-constraint calculations, (vi) block- 
ing of quasi-particles in the Equal Filling Approximation (EFA), (vii) framework for generalized 
energy density with arbitrary density-dependences, and (viii) shared memory parallelism via 
^ ■ OpenMP pragmas. 



PACS numbers: 07.05.T, 21.60.-n, 21.60.Jz 

NEW VERSION PROGRAM SUMMARY 

Title of the program: HFBTHO v2.00c 
Catalogue number: .... 

Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland (see 
application form in this issue) 

Reference in CPC for earlier version of program: M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz, 
P. Ring, Comput. Phys. Commun. 167 (2005) 43-63. 



Catalogue number of previous version: ADFL_v2_l 
^E-mail: schunckl@llnl.gov 
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Licensing provisions: GPL v3 

Does the new version supersede the previous one: Yes 

Computers on which the program has been tested: Intel Pentium-Ill, Intel Xeon, AMD- Athlon, 
AMD-Opteron, Cray XT5, Cray XE6 

Operating systems: UNIX, LINUX, Windows^P 

Programming language used: FORTRAN-95 

Memory required to execute with typical data: 200 Mwords 

No. of hits in a word: 8 

Has the code been vectorised?: Yes 

Has the code been parallelized? : Yes 

No. of lines in distributed program: 10 488 

Keywords: Hartree-Fock; Hartree-Fock-Bogolyubov; Nuclear many-body problem; Skyrme in- 
teraction; Self-consistent mean field; Density functional theory; Generalized energy density func- 
tional; Nuclear matter; Quadrupole deformation; Octupole deformation; Constrained calcula- 
tions; Potential energy surface; Pairing; Particle number projection; Nuclear radii; Quasiparticle 
spectra; Harmonic oscillator; Coulomb field; Transformed harmonic oscillator; Finite tempera- 
ture; Shared memory parallelism. 

Nature of physical problem 

The solution of self-consistent mean-field equations for weakly-bound paired nuclei requires a 
correct description of the asymptotic properties of nuclear quasiparticle wave functions. In 
the present implementation, this is achieved by using the single-particle wave functions of the 
transformed harmonic oscillator, which allows for an accurate description of deformation effects 
and pairing correlations in nuclei arbitrarily close to the particle drip lines. 

Method of solution 

The program uses the axial Transformed Harmonic Oscillator (THO) single-particle basis to 
expand quasiparticle wave functions. It iteratively diagonalizes the Hartree-Fock-Bogolyubov 
Hamiltonian based on generalized Skyrme-like energy densities and zero-range pairing interac- 
tions until a self-consistent solution is found. A previous version of the program was presented 
in: M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz, P. Ring, Comput. Phys. Commun. 167 
(2005) 43-63. 

Summary of revisions 

1. The modified Broyden method has been implemented, 

2. Optional breaking of reflection symmetry has been implemented, 

3. The calculation of all axial multipole moments up to A = 8 has been implemented. 
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4. The finite temperature formalism for the HFB method has been implemented, 

5. The linear constraint method based on the approximation of the Random Phase Approx- 
imation (RPA) matrix for multi-constraint calculations has been implemented, 

6. The blocking of quasi-particles in the Equal Filling Approximation (EFA) has been imple- 
mented, 

7. The framework for generalized energy density functionals with arbitrary density-dependence 
has been implemented, 

8. Shared memory parallelism via OpenMP pragmas has been implemented. 

Restrictions on the complexity of the problem 
Axial- and time-reversal symmetries are assumed. 

Typical running time 

Unusual features of the program 

The user must have access to (i) the LAPACK subroutines dsyevd, dsytrf and dsytri, and 
their dependencies, which compute eigenvalues and eigenfunctions of real symmetric matrices, 

(ii) the LAPACK subroutines dgetri and dgetrf, which invert arbitrary real matrices, and 

(iii) the BLAS routines dcopy, dscal, dgemm and dgemv for double-precision linear algebra 
(or provide another set of subroutines that can perform such tasks). The BLAS and LA- 
PACK subroutines can be obtained from the Netlib Repository at the University of Tennessee, 
Knoxville: http : / /netlib2 . cs . utk . edu/. 

LONG WRITE-UP 



1 Introduction 

The method to solve the Skyrme Hartree-Fock-Bogolyubov equations in the transformed har- 
monic oscillator basis was presented in [1]. The present paper is a long write-up of the new 
version of the code HFBTHO. This extended version contains a number of new capabilities 
such as the breaking of reflection symmetry, the calculation of axial multipole moments, multi- 
constraint calculations and the readjustment of the corresponding Lagrange parameters using 
the cranking approximation of the RPA matrix, the blocking prescription in odd-even and odd- 
odd nuclei, the finite-temperature formalism, and generalized Skyrme-like energy functionals. 

In addition to releasing a new version of the solver for general applications in nuclear science, 
the goal of this paper is to establish a number of precise benchmarks for nuclear structure calcu- 
lations with Skyrme functionals. To this end, we devote an entire section to comparing various 
calculations performed with the spherical HOSPHE axially- deformed HFBTHO v2.00c, and 
symmetry-unrestricted HFODD v2.49t 0, S] nuclear density functional theory (DFT) solvers. 
Also, in order to facilitate the development of future versions of HFBTHO as well as to enable 
deeper integration with future releases of HFODD, backward compatibility of input and output 
files has been broken between the version vl.66 of [1] and the current version 2.00c. Unless 
indicated otherwise, details about the methods presented in [T] still apply. 

In section ^ we review the new capabilities of the code. In section [31 we present a num- 
ber of numerical benchmarks between HFBTHO and the aforementioned DFT solvers. Such 
benchmarks are very important in view of the future development of these programs. 
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2 Modifications introduced in version 2.00c 



We present in this section the major new features added to the code between version 1.66 and 
version 2.00c. Minor improvements and bug fixes are not discussed here, the full history of 
changes can be found in the source code. 

2.1 Modified Broyden Method 

In HFBTHO v2.00c, the matrix elements of the HFB matrix are updated at each iteration 
using the modified Broyden method, instead of the traditional linear mixing of version 1.66. 
Details of the implementation, results of convergence tests, and comparisons with alternative 
implementations can be found in [3]. 

2.2 Axial multipole moments 

In HFBTHO v2.00c, the expectation value of axial multipole moments Qi = Qio = r'-Yio{6, ip) on 
the HFB ground-state is computed for all moments up to /max = 8. We recall that in spherical 
coordinates, the multipole moment Qi of order I read 

Qi{r,e,^) = r'^^^Pi{cose), (1) 

where Pi is the Legendre polynomial of order I [6] . Spherical and cylindrical coordinate systems 
are related through = p"^ + and r cos 9 = z. Recurrence relations on Legendre polynomials 
give an analytical expression for Qi{r, 2, y?) for / = 0, . . . , 8 [Ej. Multipole moments can also be 
used as constraints. In this case, the matrix elements of the Qi in the HO basis need to be 
computed. They are evaluated numerically on the Gauss-Laguerre and Gauss-Hermite nodes of 
integration used throughout the code [1]. 

2.3 Finite-temperature HFB method 

The code HFBTHO v2.00c solves the finite temperature HFB (FT-HFB) equations. The nu- 
merical implementation is similar to that of HFODD in [3]. Let us recall that the FT-HFB 
equations take the same form as the HFB equations at T = 0, only the one-body density matrix 
and pairing tensor now depend on the Fermi-Dirac occupation of quasi-particle states /i. 
Assuming axial- and time-reversal symmetry, all density matrices are real and read 

p = UfU^ + Vil - f)V^ 
K = UfV^ + V{1 - f)U^, 

with U, V the matrices of the Bogolyubov transformation. In HFBTHO, these matrices are 
block-diagonal. As in HFODD, the Fermi level A is not treated explicitly as the Lagrange 
parameter for the multipole operator Qoo alongside other multipole moments Qim- Instead, it 
is determined directly at each iteration from the conservation of particle number and is based 
on the BCS formula 

= E + i<i>^) - <W)fM ■ (3) 



4 



The BCS occupations are given by the traditional formulae 
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U, 



(4) 



{e^ — Xy + and and are the equivalent single-particle energies and 



with E^^^ 

pairing gaps, see appendix B in [7] . The Fermi-Dirac occupation factors are given by 

1 



(5) 



When using the Newton-like method to solve the equation A^(A) = N, Z for each type of particle 
at T > 0, one must now include the contribution df^/dX in the derivative of the function N{X). 



2.4 Linear constraints and the RPA method 

Multi-constraint calculations are possible in HFBTHO v2.00c. The code implements the linear 
constraint method, where the quantity to be minimized is 

E' = E-J2>^a{{QL)-QL), (6) 

a 

where Qi^ is the multipole moment operator for the constraint a and Xa is the related Lagrange 
parameter. Lagrange parameters are readjusted at each iteration according to the procedure 
presented in [8] and also used in the latest release of HFODD [3] . The philosophy of the method 
is to associate the variation of the Lagrange parameters with a first-order perturbation of the 
generalized density matrix. 

As a reminder, we start with the variations 6TZ of the generalized density matrix, which 
induce variations of the HFB matrix SH and of the Lagrange parameters 6X = {SXi, . . . ,6Xn), 
(up to first order). Neglecting the variations of the HFB matrix with respect to the generalized 
density matrix is equivalent to working at the so-called cranking approximation, and it reduces 
the HFB equation with the perturbed quantities to 

[57^,7^W]-iJ]5A47^(^Q.J=0, (7) 

a 

with TZ^^^ and T-L^^^ respectively, the unperturbed generalized density matrix and HFB Hamilto- 
nian, 6Xa the perturbation of the Lagrange parameter for the constraint a, and Qi^ the matrix 
of the constraint in the doubled s.p. basis. This equation gives the desired relation between 
STZ and 6X. The Lagrange parameter can then be readjusted at each iteration by interpreting 
the deviation 6Qi^ = {QiJ — Qi^ from the requested value Qi^ as caused by a variation of the 
generalized density matrix 671, 

6Q,^ = ^TT{Qijn). (8) 

Knowing the deviation SQi^, we obtain 671, and thereby deduce the 6Xa needed to reproduce the 
requested value. Calculations are performed in the q.p. basis, since the unperturbed generalized 
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density and HFB matrix take a very simple form. The computational cost of the method thus 
comes essentially from transforming all relevant matrices into this basis. In HFBTHO, this 
operation can be performed separately for each f2— block. The method can also be extended 
to finite-temperature in a straighforward manner by using the Fermi-Dirac occupation factors. 
Details of this extension are presented elsewhere [3] . 

2.5 Quasi-particle blocking 

Odd-even and odd-odd nuclei can now be computed by HFBTHO v2.00c using the blocking 
of quasi-particle states [9]. Because time- reversal symmetry is built into the code, the equal 
filling approximation (EFA) has to be used [10]. However, it was shown in [9] that the EFA is 
an excellent approximation to exact blocking. The identification of the blocking candidate is 
done using the same technique as in HFODD [TT] : the mean- field Hamiltonian h is diagonalized 
at each iteration and provides a set of equivalent single-particle states. Based on the Nilsson 
quantum numbers of the requested blocked level provided in the input file, the code identifies the 
index of the q.p. to be blocked by looking at the overlap between the q.p. wave-function (both 
lower and upper component separately) and the s.p. wave-function. The maximum overlap 
specifies the index of the blocked q.p. 

2.6 Generalized energy density functionals 

The kernel of the HFBTHO solver has been rewritten to enable the use of generalized Skyrme 
functionals that are not necessarily derived from an effective pseudo-potential such as the Skyrme 
force. Generalized Skyrme functionals are defined here as being the most general scalar, iso- 
scalar, time-even functional "H of the one-body local density matrix p(r) up to second-order in 
spatial derivatives of p 113]. Assuming time-reversal symmetry, such functionals thus take 
the form 

Mp] = + crWtn + Cf [p]J' + 1p]piAp, + "[p]p,V ■ J, (9) 

where t stands for the isoscalar (t = 0) or isovector (t = 1) channel, and Tt and Jt are the kinetic 
energy and spin current density in each channel. The terms C"" [p] are (possibly arbitrary) 
functions of the local isoscalar density po{r). Note that all commonly used Skyrme forces or 
functionals fall into this category because of the phenomenological density-dependent term. 
Although most Skyrme functionals have been fitted "as a force", the recent parametrizations 
UNEDFO and UNEDFl have looked at the problem more from a functional perspective [HI 
[T5j . Microscopically-derived EDF obtained, for example, from the density matrix expansion of 
effective nuclear potentials, are less trivial examples of these generalized functionals, since the 
density- dependence of the coupling constants can be significant |13] . 

In the current version, the code only implements 2'^'^-order generalized Skyrme functionals 
and it is left to the user to code more advanced functionals. The subroutine calculate_U_parameters () 
in module UNEDF provides a general template for such an implementation. Required are the form 
of the energy functional and at least its first partial derivatives with respect to the isoscalar po 
and isovector pi density matrices. Second-order partial derivatives are also necessary to compute 
nuclear matter properties. 
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2.7 Shared memory parallelism with OpenMP 

To facilitate large-scale applications of the HFBTHO solver on leadership class computers, the 
original source file has been split into a DFT solver kernel and a calling program. In version 
2.00c, we have also parallelized a number of time-intensive routines using OpenMP pragmas. 
The routine hf bdiag diagonalizes the fi— blocks of the HFB matrix: these diagonalizations are 
now done in parallel. The routine coulom computes the direct Coulomb potential Vc{ri, at 
the first iteration: this step is carried out in parallel but saves time at the first iteration only. The 
routine gamdel reconstructs the HFB matrix in configuration space for each f2— block by com- 
puting on-the-fiy the various one-dimensional integrals that define the matrix elements: shared 
memory parallelism is implemented for the outermost loop corresponding to the blocks. 
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Number of threads 



Figure 1: Time per iteration of a HFB calculation in a full spherical basis of A^max = 20 shells 
as a function of the number of threads. Tests were performed on a Intel 8-core Xeon E5-2670 
at 2.6 GHz using the Intel Fortran Compiler 12.1 and the MKL library 10.3.8. 

Figure [1] shows the performance improvement using multi-threading. We note that version 
2.00c is slightly slower than version 1.66 if only one thread is used. This additional overhead 
comes from the calculation of the densities and fields required for generalized Skyrme functionals, 
combined with the use of the Broyden method (which uses additional linear algebra at each 
iteration). Let us emphasize that HFBTHO makes use of BLAS and LAPACK routines and 
benefits from a threaded implementation of these libraries. Nested parallelism must be supported 
by the compiler. 

3 Benchmarks and Accuracy 

There exist several comparisons of HFBTHO with other DFT solvers in the literature in both 
even-even ([16l[T^) and odd nuclei ([Hill]). In some cases, these benchmarks compare different 
approaches to solving the HFB equations, in others the emphasis is put on validation of the 
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solver. Here, we want to gather in one place a comprehensive set of validation and performance 
evaluations that can be used as reference in later developments of the code. 



3.1 Benchmarks in spherical symmetry: ^"^Pb and Sn 

In spherical symmetry, HFBTHO was benchmarked against the spherical DFT solver HOSPHE 
and the symmetry- unrestricted DFT solver HFODD. We study the Hartree-Fock approximation 
in 208p]^ g^j^^ ^Yie Hartree-Fock-Bogolyubov approximation with density-dependent delta pairing 
forces in ^^°Sn. 



3.1.1 Hartree-Fock computation of Pb 

In table [H we present the results of the benchmarks between the three solvers for the spherical 
HF point in ^osp]^ ^j^g SLy5 Skyrme functional of [19]. Calculations were performed in 
-^max = 16 full spherical oscillator shells with a constant oscillator length of 6 = 2.0 fm~^. In 
HFBTHO, the number of Gauss- Legendre points for the integration of the Coulomb potential 
was iVLeg = 30, see also section 1331 below. The Skyrme energy is defined from HFBTHO outputs 
as 

^Skyrme -|- ^surface _j_ JTjSO _j_ ^tensor (10) 



HOSPHE 



HFBTHO 



HFODD 



Without Coulomb 



<i [MeV] 
4fn [MeV] 

^Skyrme [MeV] 

^so [MeV] 



I rms 

^(P) 
' rms 



[fm] 
[fm] 



2614.806852 
1438.160641 
-6498.897708 
-109.091691 
5.519846 
5.249812 



2614.806850 
1438.160644 
-6498.897712 
-109.091691 
5.519846 
5.250015 



2614.806850 
1438.160641 
-6498.897709 
-109.091691 
5.519846 
5.250015 



With Coulomb 



[MeV] 
kin [MeV] 



kin 



^Skyrme [MeV] 

^so [MeV] 
[MeV] 

7-,(cxc) 

-^Cou 

Trms [fm] 

A^s [fm] 



[MeV] 



2535.409641 
1340.663301 
-6306.660514 
-98.293331 
829.308809 
-31.312656 
5.608237 
5.448516 



2535.409574 
1340.663215 
-6306.660344 
-98.293318 
829.309073 
-31.312655 
5.608237 
5.448712 



2535.409639 
1340.663314 
-6306.660527 
-98.293331 
829.308777 
-31.312655 
5.608237 
5.448711 



Table 1: Benchmark of the three solvers HOSPHE, HFBTHO and HFODD for a spherical 
Hartree-Fock calculation in ^'^^Pb with the SLy5 Skyrme functional in a full spherical basis of 
-^max = 16 shells with oscillator length b = 2.0 fm~^. See introduction of section|3]for additional 
numerical information. 
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Without Coulomb potentials included, we note that the difference with HFODD is not greater 
than 3 eV on energies (-Eskyrmc), and the radii agree up to at least the 6*^ digit. Comparisons 
with HOSPHE show the difference in energies is less than 4 eV, while the proton radius differs 
by 0.0002 fm. It is not clear if this tiny difference comes from HOSPHE or HFBTHO. Let us 
note that the kinetic energy contribution to the total energy is probably the most sensitive to 
the details of the numerical implementation. With the Coulomb potential included (both direct 
and exchange), the discrepancy on the total energy is at most 200 eV. 

3.1.2 Hartree-Fock-Bogolyubov computation of ^^°Sn with the Lipkin-Nogami pre- 
scription 





HOSPHE 


HFBTHO 


HFODD 




Without Coulomb 


[Me VI 


1384.055496 


1384.055496 


1384.055733 


p;(p) fMeVl 


885.705972 


885.705974 


885.705876 


-Em.,,rmo iMeV] 


-3628.923682 


-3628.923686 


-3628.924198 


Eso [Me VI 


-58.837668 


-58.837668 


-58.837805 


' rms [^^^^\ 


4.678179 


4.678179 


4.678179 


vi^fi^ ffnil 

' rms [■'"'"'■■'■J 


4.455761 


4.455761 


4.455760 


Epair [MeV] 


-12.646023 


-12.646023 


-12.645709 


A^'^) [MeV] 


0.910071 


0.910071 


0.910072 


A(P) [MeV] 


0.531364 


0.531364 


0.531338 




With Coulomb 


E^l [MeV] 


1345.226426 


1345.226372 


1345.226444 


4n! [MeV] 


837.571440 


837.571383 


837.571444 


^Skyrme [McV] 


-3538.591643 


-3538.591523 


-3538.591660 


^so [MeV] 


-48.652092 


-48.652077 


-48.652100 


eS [MeV] 


367.071214 


367.071360 


367.071184 


4o^ [MeV] 


-19.140103 


-19.140102 


-19.140103 


r^ml [fm] 


4.733892 


4.733892 


4.733892 


rills [fm] 


4.585076 


4.585610 


4.585609 


^pair [MeV] 


-11.125231 


-11.125228 


-11.125232 


A(") [MeV] 


0.864875 


0.864875 


0.864875 


A(P) [MeV] 


0.481236 


0.481236 


0.481237 



Table 2: Benchmark of the three solvers HOSPHE, HFBTHO and HFODD for a spherical 
Hartree-Fock-Bogolyubov calculation in ^^°Sn with the UNEDFO functional (thus including the 
Lipkin-Nogami) prescription with a spherical basis of A^max = 16 shells with oscillator length 
b = 2.0 fm. See introduction of section |3] for additional numerical information. 

In table El we present the results of the benchmarks between the three solvers for the spherical 
HFB point in ^^°Sn for the UNEDFO Skyrme functional of [H]. Calculations were performed 
with the same basis and integration characteristics as in the previous section. The pairing 
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channel was parametrized by a density-dependent delta-pairing force with mixed volume and 
surface features, of the general type 

l^ir^ W = Vi'--^ (l - Sir - r'), (11) 

with V^""'^^ the pairing strength for neutrons (n) and protons (p), Po{r) the isoscalar local 
density, and pc the saturation density, fixed at = 0.16 fm~^. Let us recall that in the case 
of the UNEDF parameterizations, the pairing strengths should not be adjusted by the user 
since they were fitted together with the Skyrme coupling constants. Recommended values are, 
respectively, l/J"^ = -170.374 MeV and V^"^ = -199.202 MeV. Because of the zero-range of the 
pairing force, a cut-off in the q.p. space has to be introduced, and we chose Ecnt = 60 MeV in 
this example. When compatibility with HFODD is required, this cutoff is sharp, namely all q.p. 
with E > Ecut are discarded from the calculation of the density. 



3.2 Benchmarks in even-even deformed nuclei: ^^^Pu 

Next, we present the benchmark of HFBTHO in deformed even-even nuclei against HFODD. Ac- 
curate HFB calculations in deformed nuclei require the use of a suitably deformed, or stretched, 
HO basis. Such a basis is characterized by its oscillator frequencies, Ux ojy oJz in Cartesian 
coordinates, and oj± ^ in cylindrical coordinates, as well as by the total number of states 
retained. The goal of this section is to compare basis truncation schemes between HFBTHO 
and HFODD in a realistic case. 

Stretched basis in HFBTHO - It is determined by applying the general prescription 
given in |20] to the particular case of an axially-deformed prolate basis. Let us recall that the 
starting point is an ellipsoid characterized with radii R^, Ry and Rz. Introducing the spherical 
radius i?o, we have 

Rx = Rq exp 

= R^exp { \l ^/3cos h + ] }, (12) 





Rz = i?oexp { ^/^/3cos(7) } . 

The deformation of the basis is characterized equivalently by the two parameters p and q such 
that 

with the oscillator length for coordinate p. It is then assumed that 

q = ^. P = ^- (14) 
In the special case of an axially-deformed basis (/3 > 0, 7 = 0), we find 

q = e~^v^/3, p=l. (15) 
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From the volume conservation condition {b\bz = b^), Eq. lfTSll leads to 

boe-^/', b, = bp/'. (16) 



Given a "spherical" oscillator length feg and the deformation /3 of the basis, the formula f|T6|) 
uniquely defines the HO lengths of the stretched basis. 

Stretched basis in HFODD - The starting point is a general nuclear shape parametrized 
by a surface E characterized by the deformation parameters through 



R{9,^)=Ro c{a) 



A=2 /i=-A 



(17) 



where Rq = tqA^^^, c{a) is computed to ensure volume conservation, Yx^{9, ip) are the spherical 
harmonics, and the are the deformation parameters. The surface defined by f|T71) encloses 
a volume V and the radius of this "ellipsoid" along the direction ii {=x,y,z) is determined 
according to 

Frequencies of the HO along each Cartesian direction satisfy Uq = u^UyU^ with 

l^x = i^o{RxzRyz), = (^o{Rxz/Rxy) , ^^z = l^oiRxyRxzY^^ , (19) 

with the geometrical ratios R^i, = R^/ Ry. Oscillator lengths are obtained from the usual formula 



K = ^ (20) 

Discussion - It is straightforward to see that in the particular case of a prolate ellipsoid (/3 > 
0, 7 = 0, Rxy = 1/p, Rxz = Ryz = ^/q), both HFODD and HFBTHO prescriptions to choose 
the oscillator frequencies are in principle identical. In practice, however, the determination of 
the radii from Eq. (fT8|) in HFODD produce small numerical deviations compared to the analytic 
formula fll2p . This will induce systematic differences between the HO frequencies computed in 
the two codes, which will in turn alter the selection of the basis states. Figure [2] quantifies this 
statement in an extreme case. 

Figure [2] shows the numerical difference between the two codes for the total energy in ^^*^Pu 
computed for A^max = 16 and iVgtates = 500 as a function of the deformation /3 of the basis. 
The configuration chosen was obtained by putting a constraint on the quadrupole moment 
(Q20) = 150 b and hexadecapole moment {Q^o) = 30 b^; expectation values of Qeo and Qso vary 
with the deformation of the basis. For a configuration with such large deformations, a basis with 
only 500 states and A^max = 16 is not sufficient to reach convergence. In particular, important 
intruder orbitals are missing. As a result, all physical observables depend quite significantly on 
basis parameters such as its deformation or frequencies. It is therefore a good test-bench for 
numerical comparisons and is an illustration of the worst-case scenario. 

Two sets of results, with and without Coulomb potentials included, are presented. In contrast 
to the much simpler cases of section 13.11 the difference between the two codes reach up to 400 
keV, even without Coulomb terms. This discrepancy is entirely attributable to the slightly 
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Basis deformation /3 

Figure 2: Difference between HFBTHO and HFODD total energy for a very deformed configu- 
ration in ^'^"Pu (see details in text) as a function of the axial deformation of the basis (3. 

different frequencies/lengths, the impact of which is magnified by the large deformation of the 
requested configuration combined with the relatively small size of the HO basis. As an example, 
for (3 = 0.7, the oscillator lengths are b± = 2.0803 fm and b, = 2.8114 fm in HFODD, to be 
compared with b± = 2.0596 fm and bz = 2.8682 fm in HFBTHO. The Coulomb term does not 
qualitatively change this picture. Most importantly, if HO lengths are manually enforced to be 
numerically identical in the two codes, or in the case of a spherical basis, the agreement between 
the two sets of calculations without Coulomb goes back to 1 eV level as in the previous sections. 

3.3 Benchmark in deformed odd nuclei: ^^^Ba 

The new version of HFBTHO enables blocking calculations in odd-even or odd-odd nuclei. 
Since by construction HFBTHO conserves time-reversal symmetry, the blocking prescription is 
implemented in the equal filling approximation, and the time-odd fields of the Skyrme functional 
are identically zero. In [5], a detailed comparison of blocking calculations between the HFBTHO 
and HFODD solvers was presented. Calculations were done in spherical symmetry, and identical 
oscillator lengths were used. The goal of this section is to present a benchmark result for an 
odd-even nucleus in a deformed basis. As in the previous section, we do not manually enforce 
identical oscillator lengths. Instead, we use the same basis selection rules in their respective 
implementation. 

Calculations were performed in the nucleus ^^^Ba using ^^^Ba as the even-even core, with the 
SLy4 Skyrme functional, a mixed surface-volume pairing force with Vq = —300 MeV for both 
protons and neutrons, a q.p. cut-off of E^ut = 60 MeV, and a partially-filled deformed basis 
with A'^niax = 16 shells, /3 = 0.2, and Nsta^tcs = 500. In HFODD calculations, time-odd fields have 
been zeroed. Results are presented in Table 13. 3[ 

Numerical agreement is of the order of 1 keV, with maximum deviations of up to 4 keV. Such 
an agreement is in line with the results shown in the previous two sections. Since the nucleus 
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q.p. enerj 


rj [MeV] 


q.p. 


HFBTHO 


HFODD 


[512]5/2 


-1236.982297 


-1236.985082 


[633] 7/2 


-1237.317286 


-1237.316510 


[503] 7/2 


-1236.608053 


-1236.608662 


[510] 1/2 


-1236.445416 


-1236.448800 


[521] 1/2 


-1235.723785 


-1235.726439 


[523] 5/2 


-1235.445375 


-1235.449199 


[660] 1/2 


-1236.293345 


-1236.293434 


[514] 7/2 


-1235.799335 


-1235.799806 


[651]3/2 


-1235.433514 


-1235.434053 



Table 3: Results of blocking calculations in HFBTHO and HFODD in ^^^Ba in a stretched HO 
basis with /3 = 0.2 (see text for more details). 

is not as heavy as ^^"Pu and the requested configuration is much less deformed than the one 
considered in 13. 2^ basis truncation effects are mitigated, and the small discrepancy between the 
calculated HO oscillator lengths does not have as drastic an effect as in the previous section. 
Again, we note that if identical HO lengths are manually enforced, the numerical agreement is 
of the order of a few eV as shown in [9] . 

3.4 Transformed harmonic oscillator basis: ^'^Ni 

One of the characteristic features of HFBTHO is the implementation of the transformed har- 
monic oscillator (THO) basis. We recall that the THO basis functions are generated by applying 
a local scale transformation (LST) f{TZ) to the HO single-particle basis functions. The LST 
transforms every point (p, z) by 

P 

z ^ z' 

with the scale IZ = 7l{p, z) defined locally as 



7^ = 

The LST function / is chosen in such a way as to enforce the proper asymptotic conditions 
(exponential decay) for the density, according to the general procedure outlined in [21], [22]. We 
refer to [1] for the details of the implementation of the THO method in HFBTHO. 

The purpose of this section is to complete our collection of benchmarks, by comparing the 
results obtained in the THO basis produced with HFBTHO with an independent implementa- 
tion of the method written by one of us (N. Michel) and used in particular in [231 IM]- This 
program assumes spherical symmetry and has been developed independently: comparing the 
two implementations is a particularly stringent test. To do it, we used HFBTHO to generate 



P 



7^ '■ 




(21) 



(22) 
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the LST function f{7l) and its partial derivatives on a spatial mesh TZ^ with < 7^ < 40 fm by 
steps of 0.1 fm. These functionals were then read numerically by the spherical code. 

The test was carried out in the neutron-rich nucleus ^°Ni, for the SLy4 Skyrme functional 



and a pure surface pairing force characterized by V^''"^ = V^' = —519.9 MeV with a pairing 
cut-off of Ecut = 60 MeV. Both the HO basis used to generate the THO basis and the THO 
basis itself were spherical and contained A'^max = 20 full shells. The oscillator length was fixed 
at b = 2.0 fm. We present in Table [231 various quantities that are good indicators of potential 
numerical discrepancies. 



HFBTHO N. Michel's Code 



Skyrme [McV] 

^so [MeV] 
[MeV] 
<l[MeV] 
A^'^) [MeV] 
A(°) [MeV] 

Trms [fm] 



-2349.550203 
-61.591143 
1190.985716 
-58.592674 
1.915090 
-0.195796 
4.717326 



-2349.511233 
-61.590694 
1190.964057 
-58.664803 
1.916030 
-0.196557 
4.717375 




Figure 3: Radial profile of the neutron density in ^°Ni. Black plain line: results from HFBTHO, 
red dashed line: results from the spherical code. 

Overall, the agreement between the two implementations is very good. Indeed, we recall 
that the matrix elements of the Hamiltonian in the THO basis depend not only on the LST 
but also on its derivatives df /dTZ. Using a numerically generated LST in the spherical code is. 
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therefore, bound to lead to systematic deviations. We show in Figure [3] the radial profile of the 
corrected neutron density after the LST in both HFBTHO and the spherical code. The tiny 
deviations beyond r = 8 fm are the consequence of quantizing the LST in HFBTHO, and using 
this numerical function in the spherical code instead of a native LST. 

3.5 Precision of the Coulomb potential: ^^^Pu 

In HFBTHO, the direct term VQ'^^^\r) of the Coulomb potential is computed by introducing 




which allows one to write 

^cSo?W = e^^ r^/,(r). (24) 

The advantages of this method are that it removes the singularity at r = r', avoids double 
series of loops over the coordinates r and r' implied by the direct method, and is simple to 
implement. In particular, the one- dimensional integral fl24p can be efficiently computed by 
Gauss-Legendre quadrature. However, we found that this prescription, while very accurate at 
or near the ground-state of atomic nuclei, should be fine-tuned when dealing with very elongated 
shapes such as those encountered in nuclear fission. For such extreme geometries, the evaluation 
of the Coulomb potential becomes significantly dependent on the number of Gauss-Legendre 
integration points. 




Figure 4: Error in the calculation of the direct Coulomb potential as a function of the number 
of Gauss-Legendre integration points for two different geometries. Numerical details are given 
in the text. 

In Figure HI we plot the sensitivity of prescription fjM|) on the number of Gauss-Legendre 
integration points for two very different geometries along the most probable fission pathway of 
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^^°Pu: a point near the ground-state, characterized by a quadrupole moment {Q20) = 20 b, 
labeled Ci, and a point beyond the second barrier with {Q20) = 200 b and {Q30) = 36.5 b^/^, 
labeled C2. All calculations were performed with the SkM* Skyrme functional, with a density- 
dependent mixed pairing force, a cut-off of Ecut = 60 MeV in the q.p. space, and a deformed 
HO basis with iVmax = 20 shells, oscillator length b = 2.3 fm, and iVgtates = 800 states. The 
deformation of the basis was P2 = 0.2 for configuration Ci, and 0.9 for configuration C2. The 
reference direct Coulomb energy was computed with the code HFODD using the Coulomb 
Green's function method as presented in [25] . 

It is worth noting that the Coulomb energy does not converge to an asymptotic value with 
respect to the number of Legendre integration points. Instead, practitioners should look for a 
plateau condition, the position of which depends on the geometry of the self-consistent solution. 
For small to moderate deformations, the prescription suggested in [T] (iVLeg = 30) yields an 
excellent approximation to the plateau value (within 10 keV), while for very large deformations, 
the number of integration points should be increased to up to 80 in order to reach similar 
accuracy. 

4 Input Data File 

The input data file format has been entirely changed from version 1.66 to the current version 
2.00c. The number of additional features in the new version is justified to adopt a more flexible 
format for inputs. 

4.1 Sample input file 

The new format uses Fortran namelist structure. An example is included below, 
&HFBTHO_GENERAL 

number_of .shells = 10, oscillator_length = -1.0, basis_def ormation = 0.0, 

proton_number = 24, neutron_number = 26, type_of .calculation = 1 / 
&HFBTH0_1TERAT10NS 

number_iterations = 100, accuracy = l.E-5, restart_file = -1 / 
&HFBTH0_FUNCT10NAL 

functional = 'SLY4', add_initial_pairing = F, type_of _coulomb = 2 / 
&HFBTH0_PA1R1NG 

user_pairing = T, vpair_n = -300.0, vpair_p = -300.0, 

pairing_cutof f = 60.0, pairing_f eature = 0.5 / 
&HFBTH0_C0NSTRA1NTS 

lambda_values = 1, 2, 3, 4, 5, 6, 7, 8, 

lambda.active =0, 0, 0, 0, 0, 0, 0, 0, 

expectation.values = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 / 
&HFBTH0_BL0CK1NG 

proton_blocking = 0, 0, 0, 0, 0, neutron_blocking =0, 0, 0, 0, / 
&HFBTHO_PROJECTION 
switch_to_THO = 0, projection_is_on = 0, 
gauge_points = 1, delta_Z = 0, delta_N = / 
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&HFBTHO_TEMPERATURE 

set_temperature = F, temperature = 0.0 / 
&HFBTHO_DEBUG 

number_Gauss = 30, number_Laguerre = 30, number_Legendre = 50, 
compatibility_HFODD = F, number_states = 500, force_parity = T, 
print_time = / 

4.2 Description of input data 

We now define the classes of input used in version 2.00c. 
Keyword: HFBTHO_GENERAL 

• number_of _shells = 10: Tlie principal number of oscillator shells N. If the basis is spherical 
(see below), it is made of the A^'states = (A^ + 1)(A^ + 2)(A^ + 3)/6 states corresponding to full 
shells. If the basis is deformed, the code searches for the lowest A^states, with possible intruder 
contributions from up to the A^max = 90 HO shell. Default: 10. 

• oscillator_length = -1.0: The oscillator length, denoted Bq in this manuscript, corre- 
sponding to the spherical basis. It is related to the HO frequency by Bq = ^Juj^mo^. If the 
basis is deformed, the code uses the constant volume condition to define the and oscillator 
lengths such that h\ = hzh\. For negative values of the code automatically sets 6o by using 
hujQ = 1.2 X 41/^1/3. Default: -1.0. 

• basis_def ormation = 0.0: The axial deformation P2 of the basis. Only axial quadrupole 
deformations are possible. Negative values correspond to an oblate basis and are allowed. De- 
fault: 0.0. 

• proton_number = 24: Number of protons for this run. Only even values are allowed, see item 
proton_blocking under keyword HFBTHO_BLOCKING for dealing with odd-proton nuclei. Default: 
24. 

• neutron_number = 26: Number of neutrons for this run. Only even values are allowed, see 
item neutron_blocking under keyword HFBTHO_BLOCKING for dealing with odd-neutron nuclei. 
Defauh: 26. 

• type_of .calculation = 1: Defines the type of calculation to be performed for this run. 
If equal to 1, standard HFB calculations will be performed. If equal to -1, the code will do 
HFB+LN, where approximate particle-number projection is handled by the Lipkin-Nogami pre- 
scription in the seniority pairing approximation following [3]. Default: 1. 

Keyword: HFBTH0_ITERAT10NS 

• number_iterations = 100: The maximum number of iterations in the self-consistent loop. 
Default: 100. 



17 



• accuracy = 1 .E-5: Iterations are stopped when the norm of the density matrix between two 
iterations, max||p(") — p("~^)||, is lower than accuracy, or the number of iterations has exceeded 
number_iterations. Default: l.E-5. 

• restart_f ile = 1: If this switch is negative, calculations will be restarted from an existing 
solution (stored in an HBFTHO compatible binary file). If it is positive, calculations will be 
started from scratch by solving the Schrodinger equation for a Woods-Saxon potential with 
(possibly) an axial deformation defined by the value of the constraint on Q21 see below. If 
values of restart_f ile are ±1, ±2, or ±3, calculations will be unconstrained. The actual value 
will change the filename of the binary file. If values of restart_f ile are ±100, ±200, or ±300, 
calculations will be constrained. Default: 1. 

Keyword: HFBTH0_FUNCT10NAL 

• functional = 'SLY4': This is a 4-letter word that indicates the Skyrme functional to be 
used. Possible values are: 'Sill', 'SKM*', 'SKP', 'SLY4', 'SLY5', 'SLY6', 'SLY7', 'SKIS', 
'SKO', 'SKX', 'HFB9', 'UNEO', 'UNEl', 'UNE2'. Defauh: 'SLY4'. 

• add_initial_pairing = F: In restart mode (see restart_f ile ), this boolean variable de- 
cides if a small number will be added to all pairing matrix elements. This option can be useful 
to ensure that pairing correlations remain non-zero even when restarting from a nucleus where 
they have collapsed, such as a doubly-magic nucleus. Default: F. 

• type_of _coulomb = 2: Chooses how the Coulomb potential is treated. If 0, both the direct 
and exchange terms are neglected. If 1, only the direct Coulomb potential is included in the 
calculation. If 2, both the direct and exchange Coulomb potentials are included, the exchange 
term being treated in the Slater approximation. Default: 2. 

Keyword: HFBTH0_PA1R1NG 

• user_pairing = T: When this keyword is set to T, some characteristics of the pairing inter- 
action can be set by the user. It is always assumed that the pairing force reads 



Parameters that can be adjusted are the value of the pairing strength for protons and neutrons 
V^'^''' (which can be different), the cut-off in energies defining the q.p. entering the calculation 
of the densities, and the type of pairing force defined by the parameter a. When this keyword 
is set to F, a pre-defined pairing force is used for each Skyrme functional 

• vpair_n = -300.0: The value of the pairing strength (in MeV) for neutrons Vq in Eq.f l25l) . 
Default: depends on the Skyrme force. 




(25) 
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• vpair_p = -300.0: The value of the pairing strength (in MeV) for protons Vq in Eq.( l25|) . 
Default: depends on the Skyrme force. 

• pairing_cutof f = 60.0: The energy cut-off (in MeV) in q.p. space: all q.p. with energy 
lower than the cut-off are taken into account in the calculation of the densities. Default: 60.0 
MeV. 

• pairing_f eature = 0.5: The factor a in Eq.f l25l) . This parameter enables one to tune the 
properties of the pairing force: If equal to 0, the pairing force has pure volume character and 
does not depend on the isoscalar density; if set to 1, the pairing force is only active at the 
surface, since in the interior, p(r) ^ pc, if set to 0.5, the pairing force has mixed volume-surface 
characteristics. Only values between and 1 are possible. Default: 0.5. 

Keyword: HFBTHO_CONSTRAINTS 

• lainbda_values = 1, 2, 3, 4, 5, 6, 7, 8: This series of 8 integers define the multipolar- 
ity of the multipole moment constraints. It is informational only and is not meant to be changed. 

• lainbda_active = 0, 0, 0, 0, 0, 0, 0, 0: This line defines which of the multipole mo- 
ments operator Qi, for / = 1,...,8, will be used as constraints. When 0, the corresponding 
multipole is not used as constraint, when 1 it is used. Default: (/ 0, 0, 0, 0, 0, 0, 0, 
/) (unconstrained calculations). 

• expectation_values = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0: This line complements 
the preceding one by specifying the value of the constraint for each multipolarity /. Internally, 
the units for the multipole moment of order / are lO'xfm'. Example: In order to obtain a 
constraint value of = 5 b'^/^ = 5000 fm^, the third number must be set to 20.0. Default: (/ 
0, 0, 0, 0, 0, 0, 0, /). 

Keyword: HFBTHO_BLOCKING 

• proton_blocking = 0, 0, 0, 0, 0: This group of 5 integers defines the blocking config- 
uration for protons. It takes the form 2f2, tt, iV, n^, n^, where [N,nz,nr]^l'^ is the traditional 
Nilsson label. Recall that with time-reversal symmetry, states +Q and —Q are degenerate, and 
HFBTHO only considers states with positive values of Q by default: the sign of 2Q given above 
is not related to the actual value of but to the nucleus in which the blocking is performed. 
Specifically, 

• If 2f2 = 0, the entire group is disregarded (no blocking). 

• If 2^2 > 0, blocking is carried out in the nucleus with Z + 1 protons, where Z is the value 
given by flag proton_number. In practice, it means the resulting HFB solution corresponds 
to the (Z + 1, TV) nucleus. 

• If 217 < 0, blocking is carried out in the nucleus with Z — 1 protons, where Z is the value 
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given by flag proton_number. In practice, it means the HFB solution corresponds to the 
(Z — 1, A^) nucleus. 

Additionally, the user may request all blocking configurations within 2 MeV of the Fermi level 
in the even-even core to be computed. This automatization is activated by setting the parity tt 
to instead of ±1. For example, the line 1, 0, 0, 0, would compute all blocking configura- 
tions in the {Z + 1, A^) nucleus, while the line -7, -1, 3, 0, 3 would yield the configuration 
[303]7/2~ in the {Z — 1,N) nucleus. Refer to the example included with the program for a 
practical application. Default: (/ 0, 0, 0, 0, /). 

• neutron_blocking = 0, 0, 0, 0, 0: This group of 5 integers defines the blocking configu- 
ration for neutrons. It obeys the same rules as for the protons. Default: (/ 0, 0, 0, 0, /). 

Keyword: HFBTHO_PROJECTION 

• switch_to_THO = 0: This switch controls the use of the transformed harmonic oscillator basis. 
If equal to 0, then the traditional HO basis is used; if equal to -1, then the code first performs 
a calculation in the HO basis before automatically restarting the calculation in the THO basis 
after the local scale transformation has been determined; if 1, the code runs the calculation 
in the THO basis only. Note that the use of the THO option requires a large enough basis, 
typically with at least A^max = 20. Default: 0. 

• projection_is_on = 0: Particle number projection (after variation) is activated by switching 
this integer to 1. Default: 0. 

• gauge_points = 1: The implementation of particle number projection is based on the dis- 
cretization of the integration interval over the gauge angle. The number of gauge points is given 
here. Default: 1. 

• delta_Z = 0: If particle projection is on, HFB results will be projected on Z + 6Z, where Z 
is the actual number of protons in the nucleus and 6Z is specified here. Default: 

• delta_N = 0: If particle projection is on, HFB results will be projected on + 6N, where A^ 
is the actual number of neutrons in the nucleus and 6N is specified here. Default: 

Keyword: HFBTHO_TEMPERATURE 

• set_temperature = F: For finite-temperature HFB calculations, set_temperature must be 
set to T. Default: F. 

• temperature = 0.0: In finite-temperature HFB calculations, the value of the nuclear temper- 
ature is given here, in MeV. If set_temperature = F, but the nuclear temperature is positive, 
the code overwrites the fiag set_temperature. Default: 0.0. 
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Keyword: HFBTHO_DEBUG 



• number_Gauss = 30: Number of Gauss-Hermite integration points for integrations along the 
z-axis (elongation axis). Default: 30 (conserved parity), 60 (broken parity). 

• number_Laguerre = 30: Number of Gauss-Laguerre integration points for integrations along 
the perpendicular axis. Default: 30. 

• number_Legendre = 50: Number of Gauss-Legendre integration points for the calculation of 
the direct Coulomb potential, see section 3.8 of [1] and section [375] in this manuscript. Default: 
50. 

• compatibility_HFODD = F: This boolean flag enforces the same HO basis as in HFODD. In 
practice, it is only meaningful in deformed nuclei. Default: F. 

• number_states = 1000: When compatibility with HFODD conventions is enforced, this pa- 
rameter gives the total number of states in the basis. Default: Inactive. 

• f orce_parity = T: This boolean flag enforces the conservation or breaking of parity, e.g. de- 
pending on the multipolarity of the multipole moments used as constraints. Default: T. 

• print_time = 0: If 1, the time taken by some of the major routines will be printed in the 
output. Default: 0. 

5 Program HFBTHO v2.00c 

The program HFBTHO comes in the form of two files: 

• hf btho_200c . f 90 - Main file containing the self-contained HFBTHO solver. This file 
contains several Fortran modules, see below. 

• main_200c . f 90 - Calling program. 

The programming language of most of the code is Fortran 95, while legacy code is still written, in 
part or totally, in Fortran 90 and Fortran 77. The code hfbtho requires an implementation of 
the BLAS and LAPACK libraries to function correctly. Shared memory parallelism is available. 

5.1 Fortran Source Files 

The main file hf btho_200c . f 90 contains the following Fortran modules: 

• HFBTHO_VERSION: informational module only containing the change log, 

• HFBTHO_utilities: definition of integer and real number types, 

• linear _algebra: collection of various routines dealing with interpolation. 
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• UNEDF: module computing the Skyrme-like energy density and the corresponding Hartree- 
Fock fields at a given density p, 

• HFBTHO: module storing all public variables used throughout the code, 

• HFBTHO_gauss: collection of routines and functions dealing with the integration meshes. 
Contains several Fortran 77 legacy routines, 

• HFBTHO_THO: module in charge of the THO transformation. 
The rest of the routines are not stacked into a module. 

5.2 Compilation 

The program is shipped with a Makefile that is preset for a number of Fortran compilers. The 
user should choose the compiler and set the path for the BLAS and LAPACK libraries. To 
compile, type: "make" 

5.3 Code execution 

Assuming an executable named main and a Linux system, execution is started by typing 

"./main < /dev/null >& main. out " 

The program will attempt to read the file named hfbtho_NAMELIST.dat in the current di- 
rectory. The user is in charge of assuring this file is present and readable. The code will automati- 
cally generate a binary file of the form [shape] [neutron_number] _[proton_number] . [extension] 
where: 

• [shape] is one of the letters 's', 'p', 'o', which refers to spherical, prolate or oblate shape 
respectively. The choice of this letter is left to the user through the keyword restart_mode. 
This format remains for backward compatibility, 

• [neutron_number] is a 3-integer number giving the neutron number (left-padding with 
zero if necessary), 

• [proton_number] is a 3-integer number giving the proton number (left-padding with zero 
if necessary), 

• [extension] is either 'hel' (normal HO run) or 'tel' (THO run). 
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